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SUMMARY 

i 

Linear, quadratic, and cubic isoparametric hexahedral solid elements have 
| been added to the element library of NASTRAN. rhese elements are available for 
static, dynamic, buckling, and heat-transfer analyses. Because the isopara- 
; metric element matrices are generated by direct numerical integration over the 
; volume of the element, variations in material properties, temperatures, and 
stresses within the elements are represented in the computations. In order to 
; compare the accuracy of the new elements, three similar models of a slender 
! cantilever were developed, one for each element. All elements performed well. 

As expected, however, the linear element model yielded excellent results only 
| when shear behavior predominated. In contrast, the results obtained from the 
quadratic and cubic element models were excellent in both shear and bending. 

INTRODUCTION 

New aerospace vehicle concepts, such as the Space Shuttle, have added 
impetus for the continued updating of NASTRAN with the best state-of-the-art 
finite element technoligy. In response to this need, the three-dimensional 
family of linear, quadratic, and cubic isoparametric hexahedral solid elements 
were developed for and installed in NASTRAN. These three new elements signifi- 
cantly improve NASTRAN’ s capability to solve any three-dimensional solid prob- 
lem requiring static, dynamic, buckling, and/or heat-transfer analysis. 

THEORETICAL BACKGROUND 

Hexahedron solid isoparametric elements may be used to analyze any three- 
dimensional continuum composed of isotropic or anisotropic materials. Examples 
- Include thick inserts in rocket engine nozzles, thermal protection system 
. insulations, soil structure Interaction problems, and geometrically complex 
thick-walled mechanical components such as pumps, valves, etc. These solid 
elements have only three degrees of freedom at each grid point (the three dis- 
placement components), and they may be combined with all other nonaxisymmetric 
NASTRAN elements. 

The isoparametric solid elements were first presented by Irons, Ergatoudis 
and Zlenklewicz [Refs. 1 to A] . Isoparametric solid elements employing either 
eight, twenty or thirty-two grid points have been found to be suitable to solve 
most problems (Figure 1). These elements correspond to assuming a linear, 

I quadratic, and cubic variation of displacement, respectively. Clough [Ref. 5] 

I conducted an evaluation of three-dimensional solid elements and showed that the 
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pointed out that the choice of which isoparametric element is best to use 
depends on the type of problem being solved. For problems involving shear and 
bending type deformations, the higher order elements are preferred over the 
linear elements which should be used for problems in which shear stresses pre- 
dominate. It is for this reason that all three isoparametric elements have 
been incorporated into NASTRAN. 

The governing equations for isoparametric elements are based on minimum 
energy principles. The derivation of these equations assumes a displacement 
function over the element which depends on grid point displacements only. The 
governing equations are obtained by minimizing the Potential Energy which is 
evaluated in terms of these displacement functions. 


Jisplacement Functions 

The name isoparametric is derived from the fact that the interpolating or 
shape functions used to represent the deformation of the element are also used 
to represent the geometry of the element. This choice insures that the element 
displacement functions satisfy the criteria necessary for convergence of the 
finite element analysis [Ref. 4]. Referring to the curvilinear coordinates 

C) shown in Figure 1, the rectangular basic x,y,z coordinates at any point 
in the cement are obtained from the NASTRAN basic coordinates at each of the 
n grid points by: 
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(i) 


where the are shape functions which depend on the number of grid 

points used to define the element geometry. The functions are either linear 
quadratic, or cubic, and correspond to employing two, three, or four grid points 
respectively, along each edge of the element. This choice insures that there 
are no geometric gaps between grid points. Expressions for the shape functions 
may be found in Reference 6. 


The deformation of the element is represented with the identical interpo- 
lating functions used to define the geo. cry; that is: 



( 2 ) 


where u, v and w are displacements along the x, y and z basic coordinate axes. 
The displacement functions satisfy the required convergence criterion 

of adequately representing a constant strain state, and insure interelement 
compatibility along the complete element boundary [Ref. 4]. 
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Strain-Displacement Relations 


Equation (2) may be used in the well-known strain-displacement relations for 
a three-dimensional continuum [Ref. 7] to define the strain vector {e } in terms 
of the grid point displacements: 
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In order to evaluate the strain matrix [C], the derivatives of the shape 
functions with respect to x, y, and z must be calculated. Since is 

defined in terms of 5, n and it is necessary to use the relation that 
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where [J] is the Jacobian matrix. It is easily evaluated by noting that 
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where the subscripts i, 2* ... n denote the n grid points of an element. 

Stress-Strain Relations 

The stress-strain relations for a general elastic material are 

{a} - [G e ] {e - e t ) (7) 

where {o} is the 6x1 stress vector in the basic coordinate system, [G e ] is a 
6x6 symmetric elastic material matrix, and (e t ) is the 6x1 thermal strain vector. 
This thermal strain vector is defined as 


{e t > 


{o e > 


£ Nj (5,n,;) 
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( 8 ) 


where {Og} is a vector of 6 thermal expansion coefficients, and is the tem- 
perature at the ith element grid point. 

Stiffness, Mass, and Load Matrices 

The stiffness, mass, and load matrices for the Isoparametric element are derived 
by application of the Principle of Virtual Work. These element matrices, rela- 


tive to the basic coordinate system, are given by 
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[where [K ee ] is the element stiffness matrix in the basic coordinate system, 
j [M ee ] is the mass matrix, {P^} is the thermal load vector, and {PgP} is the 
j pressure load vector derived from surface pressures on each of the six faces 
i of the solid element. |j| is the determinant of the Jacobian matrix, and [N] 
is a matrix of the isoparametric shape function defined by 
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P_£ is the uniform normal pressure (positive outward) applied to the face of 
the element where £ * -1; P_^ is the pressure applied to the face where 

n * -1. etc.; and {J^}, { J~^} , and {J~^} are the first, second, and third col- 
umns, respectively, of the inverse of the Jacobian matrix. Products like 
| J | {JgM in the expression for pressure load are equivalent to a vector of 

direction cosines multiplied by a surface area scaling factor relating the curvi- 
linear coordinates to the basic coordinate system. 

The integrals in equations (9) to (12) are evaluated numerically by using the 
method of Gaussian Quadrature [Ref. 8]. In the above equations, therefore, [C], 

( J | , and [N] must be evaluated at each interior point used for numerical inte- 
gration. [G e ], {<x e }, and p can also be evaluated at each integration point. 

Thus, variations in these quantities are allowed because of, say, temperature- 
dependent material. 
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The computations for the 1st parametric elements are carried out in the 
basic coordinate system. If the global coordinate system at any grid point is 
different from the basic system, the final matrices and vectors are transformed 
into that global system. 

Stress Recovery 

The equation for calculating element stresses at any interior point of an 
isoparametric element may be obtained by combining equations (3), (7), and (8) e 
follows : 


(a) - [GJ [C] {u J - (aj 
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[C] and are functions of the element curvilinear coordinates £,r|,5 evaluated 
at the point at which stresses are desired. 

Differential Stiffness Matrix 

The differential stiffness matrix for the isoparametric solid element is 
derived by adding the energy of an initial stress state to the potential energy 
function. This additional energy is derived in Reference 9 and is given by 
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These rotations may be expressed in terms of the grid point displacements by 
using equation (1) : 
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> Substituting equation (17) into equation (15) and adding this function to the 
potential energy expression yields the differential stiffness matrix: 
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vs with the structural stiffness matrix, this integral is evaluated using the 
nethod of Gaussian Quadrature. The differential stiffness is computed in the 
>asic coordinate system and then transformed, as required, to the NASTRAN 
global system. 


IMPLEMENTATION 

Many existing NASTRAN functional modules and subroutines were modified to 
Implement the Isoparametric solid elements. Several new subroutines were also 
^dded. These modules and brief descriptions of the changes to each are listed 
^.n Table 1. The detailed description of these changes presented in Reference 10 
:an be used to augment the NASTRAN Programmer's Manual instructions, Section 6.8, 
to assist in the installation of other new elements of similar complexity. Many 
if the changes are those normally required when implementing new elements. How- 
ever, in this case, changes were alBO required in the PL0T module (to plot 
three-dimensional elements), the GP3 module (to process a new external pressure 
Load), .nd in various other modules to accomodate the large space requirements 
jf the 32- grid-point cubic element. 

It should be noted that the isoparametric elements were installed in func- 
tional modules SMA1, SMA2, and DSMG1 on an interim basis only. The element 
natrlx subroutines were designed specifically for the new Element Matrix Gener- 
ator module and will be made available with Level 16 of NASTRAN. 


EVALUATION 


I A slender cantilever beam model was chosen to evaluate the performance of 
he three new isoparametric solid elements in NASTRAN. This model was chosen 
or two reasons: (1) theoretical solutions are well known, and (2) solid finite 
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elements characteristically do not perform well when used to model structures 
which exhibit predominant bending behavior. 

Three models were prepared as shown in Figures 2 3, and 4, one with each 

of the three elements: IHEX1, IHEX2, and IHEX3, the linear, quadratic, and 

cubic elements, respectively. All three beam models had a length (L) of 3.66m 
(144 in.) and a uniform rectangular cross section with depth (D) of 0.61m 
(24 in.) and width (W) of 0.30m (12 in.). The same uniform material properties 
shown in Table 2 were used for all three models. Static, normal modes, and 
buckling analyses were performed for each of the beam models. 


Statics 

For the static analyses, all de^-ees of freedom at the base of the beam 
(z ■ 0) were completely fixed. All three models were subjected to the same 
four loading conditions: 

1. Linear thermal gradient (Y-direction) 

T - 322.04 K at Y - 0 (120° F at Y - 0) 

T - 188.71 K at Y - 0.61m (-120° F at Y - 24 In.) 

2. Uniform temperature rise 

AT - 55.56 K (100° F) 

3. Compressive axial pressure (Z-direction) 

P z - -2.954 x 10* N/m 2 at Z » 3.66m (-42837 psi at Z - 144 in.) 

4. Transverse pressure (Y-direction) 

P Y - 6.895 x 10 s N/m 2 at Y - 0 (100 psi at Y - 0) 

The results for the tip displacements are summarized in Table 3, where the 
computed solutions are compared with the theoretical solutions. The maximum 
error for the linear IHEX1 element was 10. 3% for the transverse pressure load. 
For the quadratic and cubic elements, IHEX2 and IHEX3, the maximum errors of 
4.5% and 3.5%, respectively, occurred in the solutions for the thermal gradient 
load. For the transverse pressure load, the errors were 1.6% for the IHEX2 
element model, and 1.1Z for IHEX3. Thus, the higher order isoparametric sclid 
elements perform very well when used to model the bending behavior of this beam. 

Normal Modes 

In the normal mode analyses, the same single point constraints were 
applied to all three models in the following manner: All Z components of dis- 

placement in the plane Z ■ 0 and all Y components along the line Z • 0, 

Y ■ 0.30m (12 In.), were fixed. For the IHEXi and IHEX3 models only, all X 
components along the line Z • 0, X • 0, were fixed. For the IHEX2 model only, 
the X components along the line Z ■ 0, X - -0.15m (6 in.), were fixed. This 
system of constraints was chosen to allow dilatation at the base of the beam. 

The particular set of constraints used for the IHEX2 model has the additional 
advantage of symmetry. 
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The inverse power method was used to extract the first three normal modes 
of each model. The results for the natural frequencies are summarized In 
Table 4. The computed natural frequencies for the IHEX2 and IHEX3 ir Jels are 
within 3.0 per cent of the theoretical solution. The natural frequency for the 
IHEX1 model is 2.7% off for bending in the Y-direction, but it is off by 18.3% 
and 15.8% for the two bending modes about the X-dircction. These er: s are 

probably caused by an insufficient nunber of elements through the width of the 
beam in the X-direction. Using a smaller mesh size with more IHEX1 elements 
would improve these results at the expense of increased computer costs. This 
problem, therefore, serves to demonstrate even more clearly the superiority of 
the IHEX2 and IHEX3 elements over the IHEXI element for modeling the bending 
behavior of structures. 

All the computed mode shapes for all three models showed excellent corre- 
lation with the theoretical solution [Ref. 11]. Comparative plots of the mode 
shapes are not Included in this paper because there would be no visible dis- 
tinction between computed and theoretical solutions. 


Buckling 

Each of the three beam models was used to compute the critical buckling 
load for axial pressure. The same system of constraints used to compute normal 
modes was used to compute the axial pressure buckling load. The applied pres- 
sure on the end of the beam was -2.954 x 10* N/m 2 (42,837 psi). This amounts 
to a tctal applied force of -5.406 x 10 7 N (-1.234 x 10 7 lb), which is equal to 
the theoretical critical load for buckling in the X-direction. Therefote, the 
fundamental eigenvalue for buckling should have been unity. 

Again, the inverse power method was used to extract the three lowest buck- 
ling modes. The results for the buckling eigenvalues X are presented In 
Table 5. The IHEX2 and IHEX3 elem»**t results are excellent. They are within 
0.7% of the theoretical solution. The eigenvalue for the IHEXI element model 
is in error by lens than 10% for buckling in the Y-direction. However, it is 
off by more than 40% for both buckling modes in the X-direction. This situation 
is siuilcr to that of the normal mode problem for the IHEXI element model. 

Again, it is probably due to the lack of an adequate number of elements through 
the width of the beam in the X-direction. 

as was the case for thn normal modes problem, the mode shapes computed by 
NASTRAN for buckling were very close to the theoretical shapes. Thus, no plots 
comparing computed shapes with theoretical shapes are Included in this paper. 

CONCLUDING REMARKS 

All three isoparametric solid elements produced good results for static, 
normal mode, and buckling analyses. As expected, the linear element results 
shoved that it is best used when sheer behavior predominates. The superiority 
of the quadratic and cubic elements was confirmed by the excellent resul'.s 
obtained in both the bending and the shear behavior of a cantilever bean model. 
Therefore, the implementation of these three isoparametric solid elemer's, which 
provide for variations in both material properties and stresses throughout the 
element, does greatly enhance the total modeling capability of NASTRAN. 
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TABLE 1. FUNCTIONAL MODULE MODIFICATIONS TO IMPLEMENT 
ISOPARAMETRIC SOLID ELEMENTS 


IFP 

GP2 

PLTSET 

PL0T 

GP3 


New Bulk Data cards were added 


Array sizes were increased to accommodate elements 
with 32 grid points 


Array sizes were increased to accommodate elements 
with 32 grid points 

Capability for plotting solid elements was 
implemented 


Processing of the isoparametric element pressure 
card was implemented 


TA1 - Capability to append grid point temperatures to 
EST/ECPT entries was implemented 


SMA1 - Stiffness and conductance matrix generation for 
the new elements was implemented 


SMA2 - Mass and capacitance matrix generation for the new 
elements was implemented 


SSG1 - Load vector generation for thermal and pressure 
loads on the new elements was implemented 


DSMG1 - Differential stiffness matrix generation for the 
new elements was implemented 


SDR2 - Stress calculations for individual grid points of 
the new elements was implemented 

0FP - Stress printout formats for the new elements were 
Implemented 
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TABLE 2. MATERIAL PROPERTIES OF THE CANTILEVER BEAM MODELS 


Description 


Value (SI) 


Value (English) 


E Young's modulus 

V Poisson's ratio 


2.068 x 10 11 N/m 2 I 30 x 10 b lb/in" 


a Coef. of thermal expansion 2.570 x 10* 5 — m yy 1.428 x 10” 5 rr~ 

m - k xn - 


p Mass density 


20.86 kg/m 3 


in - F 

7.535 x 10" 4 lb/in 3 


TABLE 3. COMPARISON OF TIP DEFLECTIONS FOR NASTRAN AND THEORETICAL 
SOLUTIONS FOR FOUR STATIC LOADING CONDITIONS 


Description 


Thermal Gradient 


NASTRAN Solutions 

Theoretical 

IHEX1 1 

Model 

IHEX2 Model 

IHEX3 Model 

Solution* 

Defl, , 
cm 

Error, 

% 

Defl. , Error, 
cm % 

Defl. , Error, 
cm % 

Tip Defl., 
cm 

3.668 

2.5 

3.932 4.5 

3.894 3.5 

3.762 

.5367 

2.8 

.5344 2.3 

.5304 1.6 

.5222 

-.5179 

0.8 

-.5187 0.7 

-.5199 0.4 

-.5222 

.3612 

10.3 

.3965 1.6 

.3985 1.1 

.4028 


♦Theoretical Solutions 
Load Case 1 Load Case 2 

. frvr 2 


Y 2D 


6 7 - aATL 


Load Case 3 


Z E 


Load Case 4 


A . 3P * L 

* 3 1 + 2 

2ED J 5L Z 


6y ■ 3.762 cm 6y " *5222 cm 


6y *»-.5222 cm 


.4028 cm 


m 





TABLE 4. COMPARISON OF NATURAL FREQUENCIES FOR NASTRAN AND 

THEORETICAL SOLUTIONS 





NASTRAN Solutions 


Theoretical 

Mode 

Description 

IHEX1 

Model 

IHEX2 

Model 

IHEX3 

Model 

Solution, 

Freq . , 
cps 

Error, 

% 

Freq. , 
cps 

Error, 

% 

Freq. , 
cps 

Error , 

% 

cps 

[Ref. 11] 

i 

First Bending Mode 
in the X-Direction 

22.0 

18.3 

18.6 

0 

18.6 

0 

18.6 

2 

First Be ;ding Mode 
in the Y-Direction 

38.3 

2.7 

36.5 

2.1 

36.5 

2.1 

37.3 

3 

Second Bending 
Mode in the 
X-Direction 

135.3 

15.8 

114.3 

2.1 

113.3 

3.0 

116.8 


TABLE 5. COMPARISON OF BUCKLING EIGENVALUES FOR NASTRAN AND 

THEORETICAL SOLUTIONS 




NASTRAN Solutions 

Theoretical 

Mode 

Description 

IHEX1 Model 

IHEX2 

Model 

IHEX3 Model 

Solution 

X 

Error, 

X 

X 

Error , 
% 

X 

Error , 
% 

A 

[Ref. 12] 

1 

X-Direction 

1.406 

40.6 

1.002 

.2 

1.001 

.1 

1.0 

2 

Y-Direction 

4.391 

9.8 

3.981 

.5 

3.979 

.5 

4.0 

3 

X-Direction 

12.809 

42.3 

9.037 

.4 

8.934 

.7 

9.0 
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FIGURE 1. THREE ISOPARAMETRIC ELEMENTS 






3.66 m — 


r 0.30 m 


FIGURE 2. IHEX1 MODEL — 216 ELEMENTS AND 364 GRID POINTS 
MATRIX ORDER (g-SET) = 1092, SEMI -BANDWIDTH * 102. 


FIGURE 3. IHEX2 MODEL -- 36 ELEMENTS AND 275 GRID POINTS 
MATRIX ORDER (g-SET) = 825, SEMI-BANDWIDTH - 156. 


FIGURE 4. IHEX3 MODEL ~ 8 ELEMENTS AND 148 GRID POINTS 
MATRIX ORDER (g-SET) - 444, SEMI-BANDWIDTH - 132 


